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We present a pedagogical introduction to the theoretical framework of 
the Color Glass Condensate (CGC) and the McLerran-Venugopalan (MV) 
model. We discuss the application of the MV model to describe the early- 
time dynamics of the relativistic heavy-ion collision. Without longitudi- 
nal fluctuations the classical time evolution maintains boost invariance, 
while an instability develops once fluctuations that break boost invariance 
are included. We show that this "Glasma instability" enhances rapidity- 
dependent variations as long as self-interactions among unstable modes stay 
weak and the system resides in the linear regime. Eventually the amplitude 
of unstable modes becomes so large that the growth of instability gets sat- 
urated. In this non-linear regime the numerical simulations of the Glasma 
lead to turbulent energy flow from low-frequency modes to higher-frequency 
modes, which results in a characteristic power-law spectrum. The power 
found in numerical simulation of the expanding Glasma system turns out 
to be consistent with Kolmogorov's —5/3 scaling. 



1. Introduction 

Relativistic heavy-ion collision experiments have aimed to create a new 
state of matter out of color-deconfined particles, i.e. a quark-gluon plasma 
(QGP) in extreme environments in the laboratory. Presumably it was only 
up to ~ 10^^ s after the Big Bang when the Early Universe was still hot 
enough to realize the QGP in nature. Experimental data on intrinsic prop- 
erties of the QGP suggest that this new state of QCD matter found in the 
heavy-ion collision is not a weakly-coupled plasma but rather a strongly- 
coupled fluid. The hydrodynamic description of the time evolution has 
successfully reproduced the measured particle distributions, in particular, 
the azimuthal distribution of emitted particles in non-central collisions. The 
great success of the hydrodynamic model is a strong evidence for thermal- 
ization. After thermal equilibrium is achieved, the time evolution of the 
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QGP is somehow under theoretical control, though some uncertainties re- 
main in the determination of the equation of state and the implementation 
of dissipative effects. 

In contrast to the hydrodynamic regime after thermalization, our under- 
standing is still quite limited about the early-time dynamics toward thermal- 
ization. This kind of problem is in general one of the most difficult physics 
challenges. One would naturally anticipate that the system may form a 
turbulent fluid during transient stages right after the collision. In fact, tur- 
bulence is a common phenomenon in our daily life whenever the Reynolds 
number exceeds a certain threshold. Nevertheless, it poses a very difficult 
theory problem even today. It is widely known that a renowned physicist, 
Richard Feynman described turbulence as "the last great unsolved problem 
of classical physics" and it is indeed so also in the context of the QGP study. 

A good news is that we already have a powerful theoretical tool to inves- 
tigate the very early-time dynamics of the high-energy heavy-ion collision. 
One may well think at a first glance that any microscopic description of 
nucleus-nucleus collisions is too complicated to handle in terms of the QCD 
first principle. This is truly so unless the collisional energy is sufficiently 
high. Actually microscopic information of nucleus before collision is far 
from simple on its own. The high-energy limit, however, allows for drastic 
simplification that makes the calculation feasible. In the Regge limit, pre- 
cisely speaking, the strong interactions exhibit totally different characteris- 
tics from low-energy hadron physics. In this particular case the cm. energy 
scale s is infinitely larger than other energy scales such as the transferred 
momentum squared t. Then, although the strong coupling constant is 
small due to asymptotic freedom with large s, a resummation is required for 
a series of non-small terms (aglnx"^)" where Bjorken's x ~ t/s. In terms 
of the Feynman diagram this resummation represents quantum processes to 
emit softer gluons successively. Intuitively, as one goes to larger s and thus 
goes to x — 5x down from x, one should consider more radiated gluons within 
a bin of x to x — 6x. Naturally the wave function of nucleus is x-dependent 
and we should expect more and more gluons inside at smaller x. 

Eventually it should be the most suitable to treat gluons as coherent 
fields rather than particles once the gluon density is high enough. This is 
reminiscent of photons in the Weizsacker- Williams approximation. In the 
first approximation in the high-energy limit of QCD, therefore, the clas- 
sical fields are appropriate ingredients in theoretical computations, which 
is as a consequence of quantum radiations. Such a classical description of 
high-energy QCD is called the Color Glass Condensate (CGC) [I]. In this 
way, the initial condition for the relativistic heavy-ion collision should be 
formulated by means of the CGC theory. As long as the gluon distribution 
function stays large, the CGC picture holds, and it is finally superseded 
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by a particle picture of plasma. Some unnamed physical state between the 
CGC and the QGP was given a name of Glasma as a mixture of "glass" 
and "plasma" [2j. The Glasma time evolution thus follows the classical 
equations of motion and hopefully it should be transformed smoothly into 
a hydrodynamic regime. In this sense the Glasma should play a central role 
to figure out what the initial condition for the hydrodynamic model should 
look like. 

Because QCD or the Yang- Mills theory involves gluonic self-interactions, 
it is generally hard to find an analytical solution of the classical equations 
of motion except for some simple situations. Hence, one has to resort to 
numerical methods to go into quantitative estimates, and the numerical im- 
plementation has been established. In view of numerical outputs, however, 
there is no indication seen in the Glasma evolution toward thermalization. 
If all quantum fluctuations are completely frozen in particular, the classical 
dynamics respects boost invariance. The coordinate rapidity r] is simply 
shifted under a boost on the system, that means that the boost-invariant 
system is insensitive to r/. The pressure is then highly anisotropic depending 
on the beam-axis direction or on the transverse plane. This makes a sharp 
contrast to thermalized matter in which the pressure is isotropic. 

Quantum (or structural - see Sec.[S]) fluctuations including ?7-dependence 
break boost invariance explicitly. Interestingly enough, it was discovered 
that a small modulation along the longitudinal direction grows up expo- 
nentially as a function of time and r^-dependent modes show instability 
behavior, which is sometimes referred to as the "Glasma instability" [3j. It 
is of paramount importance to understand the nature of the Glasma insta- 
bility to fill in a missing link between the CGC initial state and the initial 
condition for hydrodynamics. There are several theoretical attempts to ac- 
count for qualitative features of the Glasma instability [H El El [3 [H [9]. 
Here, instead of doing so, we will think about subsequent phenomena at 
later stages; turbulence may be formed by instability growth, and then it 
is sensible to anticipate the decay of turbulence and the associated scaling 
law in the energy spectrum |10j . 

This article contains lectures on the basic facts of the MV model for 
those who are not necessarily familiar with QCD at high energy. In Sec. [2] a 
stationary-point approximation on the functional integration is introduced, 
which leads to a classical treatment of the high-energy QCD problems. The 
classical equations of motion in the pure Yang-Mills theory are further dis- 
cussed in Sec. [3l Then, some numerical results for the r/-independent case 
are presented next in Sec. |H and those for the case with ?7-dependent fluc- 
tuations in Sec. El Some evidence is presented for the realization of the 
power-law scaling in the energy spectrum in the non-linear regime where 
the instability stops. Section El is devoted to outlooks. 
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2. Scattering amplitude and the Eikonal approximation 

We will see the essence of the Eikonal approximation and the scattering 
amplitudes of our interest in QCD physics. We will first consider the case of 
light projectile and dense target, and proceed to the case of dense projectile 
and dense target next. 



2.1. Eikonal approximation 

Before addressing QCD application, we shall first consider a scatter- 
ing problem in non-rclativistic Quantum Mechanics. If we want to solve a 
problem of potential scattering, we should treat the Schrodinger equation. 



with the following boundary condition. 



(1) 



(2) 



at large distance (r — t- oo). The term involving /(Jl) represents the scattered 
wave and one can obtain the scattering amplitude from The lowest- 

order estimate immediately gives an expression in the Born approximation. 
In the high-energy limit, however, the scattering angle is small and the 
incident and the scattered waves interfere strongly there. In this situation 
the following Ansatz is more convenient, 



J-kz 



(3) 



which is called the Eikonal approximation in analogy with the terminology 
in Optics. The wave length of incident wave, A = 27r//s, is shorter than the 
potential range when k is large enough, and the differential equation for 
^(r) is 

r 

vPz- 7, 1- ^{f) 



2m 



tjj{r) = 



(4) 



with V = hk/m. Then, in the high-energy limit, the second term is negligible 
as compared to the first term, which is easily integrated out to give. 



■0(r) = exp 



-- r dz'V{x,y,z') 

V J-oo 



(5) 



In the gauge theory dz'V{z') is replaced by dz'A^j^z'). Let us consider the 
scattering of the target particles moving at the speed of light in the positive- 
z direction and the projectile particles in the negative-z direction. Then, in 
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Fig. 1. Schematic figure to show scattering between the target partons along the 
^'''-axis and the projectile partons along the a;~-axis. 



general, in the Eikonal approximation, the scattering matrix is 



{pt} {pt} {pp} {pp} 



with the Wilson lines corresponding to the Eikonal phase ([5]), 



W{x±,x ) = V exp 
V{x±,x~^) = V exp 



(7) 



in the light-cone coordinates; = (z + t)/\/2 and x^ = {z — t)/\/2, as 
sketched in Fig. [TJ The momenta conjugate to x"*" and x~ are the energy 
p~ = {E — p^)/\/^ and the longitudinal momentum = [E +p^)/-v/2- In 
Eq. ([6]) the weights, Wa;[/3t] and W2,.'[/0p], represent the wave functions of 
the target and the projectile at the Bjorken variables x and x', respectively. 
Since the hard particles (p+ > xP^ with being the total longitudinal 
momentum of the target) are included in the wave function, the functional 
integration in (• • • ) should contain the softer gauge fields with p+ < xP^ . 
In the following subsections let us discuss how to make an approximation 
on this scattering amplitude. 
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2.2. Light projectile and dense target 

For the simplest example let us take the projectile as a color-dipole, i.e. 
the scattering amplitude is then, 

'S'dipolc ~ ' 

= Y,y^x[pt] ! P.4y(a;x)l/^(yx)e'-^YM[A]+i5eou.cc[pt,Ty] _ (g) 
The product of W's can be re-expressed on the exponential [HI [12] as 

5source[pt, ^1 = ^ / [p^ InW] ^ - J d^X p^A^ , (9) 

where the last expressions is an approximation valid for sufficiently large pf 
Then, with large pt, the functional integration in Eq. ([8]) can be estimated 
by means of the stationary approximation at the solution of 



SSym 



= S^~pt . (10) 

A=A 



The solution of the above classical equations of motion thus represents the 
contribution from soft gluons with p+ < xP~^. Then, finally, the dipole 
scattering amplitude is 



^dipoic ^ Yl V{x^)V\y^] 

{pt} 



(11) 

A=A[pt] 



This expression is easily generalized to an arbitrary operator 0[A] and 

mA])),,^Y.^MO[A[pt]]. (12) 

{pt} 

Once the x dependence in Wa;[/Ot] is known [13l [ll], small-x evolution is 
deduced for {{0[A]))p^ in general and in this way the BFKL equation up to 
the quadratic-order of pt, and besides, the BK [IS] and JIMWLK equations 
including full-order of pt are derived. 
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2.3. Dense projectile and dense target 

The discussions so far are quite generic but the above stationary-point 
approximation needs slight modifications when not only the target parton 
density pt but also the projectile pp is large as in the situation of the rela- 
tivistic heavy-ion collision. Then, the scattering amplitude reads 

^dense-dense = ^ ^M^APp] IvAe'^^^^^^^+'^—^P^'^^'P^'''^ (13) 
{Pt,Pp} 



with the source action given approximately as 



5sourcebt, W; Pp, V] - j d^x {p^A- + p^A+) . (14) 

This time, the stationary-point is shifted also by the effect of the presence 
of pp and it is determined by the classical equations of motion with two 
sources, 

^^"""^ S^-pt + S^^+pl. (15) 



A=A 



Using the solution A of the above classical equations of motion, one can 
obtain the general formula similar to previous Eq. (jl2p as 



((o[^]»pt,pp ^ E mptmApp]o[A[puPp]]. (16) 

{pt,pp} 

In what follows we discuss how to solve this above Eq. to investigate 
the early-time dynamics of the heavy-ion collision. 

Then, there are two ingredients necessary to estimate physical observ- 
ables using Eq. (fT6]l . One is the solution of Eq. (fT5]) , which is impossible to 
find in an analytical way unfortunately. The other is to figure out the wave 
functions Wx [pt] and Wx' [pp] , which is again impossible to do so by solving 
QCD exactly. Given an initial condition at a certain x, in principle, the 
evolution equation such as the JIMWLK equation leads to the wave func- 
tion at any x. The theoretical framework of such a description of scattering 
processes with non-linearity of abundant gluons is called the Color Glass 
Condensate (CGC). The CGC theory is not a phenomenological model but 
an extension of conventional perturbative QCD with resummation in a form 
of background fields. In any case, however, the initial condition for small-x 
evolution is necessary for actual computations. This part needs an Ansatz 
as explained in the next section. 
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3. Equations of motion and the MV model 

Here we will see that the analytical solution is written down for one- 
source problem. Although it is impossible to give an analytical formula for 
the solution of two-source problem, the initial condition on the light-cone 
can be specified. Also we will introduce a Gaussian approximation for the 
wave functions that defines the MV model. 



3.1. One- source problem 

Let us first consider how to solve the equations of motion (jlOp for the 
light-dense scattering. It is actually easy to find a special solution in the 
same way as the classical electromagnetism. The important point is that the 
source p\,{x±,x~) is independent of x"*" because of time dilation of particles 
moving at the speed of lighlQ. Therefore, with an assumption of 7^ 
and A~ = = 0, the problem is reduced to an Abelian one and Eq. ([TO]) 
amounts to the standard Poisson equation. The solution of the static po- 
tential therefore reads, 



A-" 



1 

9 



2 Pt(a;_L,x" 



(17) 



In later discussions it is more convenient to adopt the light-cone gauge, 
A'^ = 0, which is achieved by a gauge rotation by that solves [A'^ — 
{\g)-^d+]V^ = 0. That is, 



V\x^,x ) = Vexp 
= V exp 



15 



-15 



dz A'^{x_\_,z 



dz —^ptix±,z' 
(J I 



(18) 



After the gauge transformation by V\ only the transverse components are 
non-vanishing, 

of = A = -lgy^^V^ ' (19) 

and A'^ = A~ = 0. In the same way, for the projectile moving in the 
opposite direction to the target, the equations of motion have a solution as 



(p) 



1 

iff 



WdiW^ 



(20) 



^ Even though the time dependence is frozen at the speed of hght, non-commutativity 
of color charges needs distinction by the order of interactions along the a;+-axis, which 
introduces a label that plays the role of time. Such "a;+-dependence" is dropped off 
in the large A'c limit only. 
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with 



W^x^.x-^) = Pexp 



-15 / d2;+--2-pp(a;_L,z+) 



(21) 



in A~ = gauge. In this manner the one-source problem is readily solvable. 
It is, however, impossible to solve the two-source problem in Eq. (jlSp as 
simply as above. 

3.2. Two- source problem 

In the presence of two sources it is most suitable to make use of the 
Bjorken coordinates that represent an expanding system. The time and the 
longitudinal variables are replaced by the proper time and the space-time 
rapidity, respectively, as 

r = \'2x+x- = _ ^2 ^ ri = -\n— = -hi(^-^\. (22) 

2 X 2 \t — z J 

The temporal gauge in this coordinate, Ar = 0, has a close connection 
to the light-cone gauge discussed in the previous subsection because Ar = 
x~A'^ + x^A~ = leads to a condition A'^ = on x"*" = and A~ = on 
x^ = 0. Therefore, the experience in the one-source problem turns out to 
be useful here. 

The equations of motion can be expressed in the Bjorken coordinates as 

I (23) 
and the conjugate momenta are defined as 

£^ = TdrAi , = ^drAr, ■ (24) 

These equations ()23p and ()24p determine the time evolution uniquely once 
the initial condition at an initial time tq is specified. 

The solutions of the one-source problems are consistent with the gauge- 
fixing condition Ar = since we found A'^ = A~ = for the target and the 
projectile both. Then af^ solves Eq. ()23p too if there is no interference from 
the projectile source. This means that af^ should be a solution in the region 

outside of the light-cone, (x+ < 0, x^ > 0), and a^^^ in {x~^ > 0, x~ < 0) 
(see Fig. [5] for illustration). 
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Fig. 2. Solution of the equations of motion (jTS]) outside of the forward hght-cone. 



Prom this fact, it is naturally understood that the initial condition at 
r = or x+ = x~ = may be a superposition of these two solutions, 

Ai = + af^ , A^ = 0. (25) 

Here, though this is a very simple Ansatz by a superposition, the conse- 
quence is quite non-trivial. The field strength associated with these gauge 
fields is 

B^ = 0, B^ = Fi2 = -ig[[4\a^i^] + [a^^\af]) , (26) 

and thus, the longitudinal component of the chromo-magnetic fields appears 
from the non-Abelian interactions. By solving the equations of motion, one 
can also find similar expressions for the chromo-electric fields [16] 

f ^ = , E"^ 



■.g{yi\af^] + [4\af]). (27) 



These field strengths stand for characteristic properties of the initial con- 
dition of the relativistic heavy- ion collisions in the CGC or the so-called 
Glasma picture, an intuitive illustration for which is displayed in Fig. O It 
is important to point out that the initial conditions at r = are indepen- 
dent of r/, namely, boost invariant. Because there is no explicit r/ in the 
equations of motion, boost invariance is kept during the time evolution. 

3.3. McLerran-Venugopalan model 

One of the simplest and reasonable Ansatze for the wave function is the 
Gaussian approximation, that is given by 



Wx[p] = exp 



(28) 
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Fig. 3. Intuitive picture of the initial condition for the heavy-ion colhsion charac- 
terized by the CGC initial conditions ([25]) and (|24l) . 

which defines the McLerran-Venugopalan (MV) model. Here z represents 
either x"*" or x~ . In this model setup characterizes the typical energy 
scale. Indeed, once the parton saturation manifests itself, any details in the 
structure are lost and only the transverse parton density ~ should be 
a relevant scale. In principle ^{z) in the MV model is to be interpreted 
as Qs in the parton saturation. With the Gaussian wave function ()28p . 
the expectation value is obtained by a decomposition into the two-point 
function that is read as 



In other words, the Gaussian approximation (|28p assumes no correlation at 
all between spatially distinct sites. Evaluating the Gaussian average with 
various functionals of p°'{x^,z) is an interesting mathematical excercise [17j . 
Especially it is feasible to evaluate the initial energy density, 



at r = using the initial conditions ()25p and ()27p and the Gaussian 
weight (j28p . The results are a bit disappointing because it involves both 
the UV and the IR divergences, which can be regularized by A (momentum 
cutoff) and m (gluon mass). Then, the initial energy density is found to be 



after some calculations [U [T8l [T9] . In the numerical simulation with dis- 
cretized grid in a finite-volume box |20l [2T| there are natural UV and IR 
cutoffs incorporated from the beginning. The IR cutoff is, however, not 
included as a gluon mass but originates from the system size L. 




(30) 




(31) 
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Fig. 4. Left: An example of the initial distribution of A(a;^) as a solution of the 
Poisson equation with random color source. Right: The corresponding gauge field 
a^ = -(l/iff)y5.yt. 

4. Numerical method and the boost-invariant results 

The model parameters should be fixed first. In the case of the gold-gold 
collision at y/s^ = 200 GeV, the empirical choice is g = 2 and g^^L = 120, 
which correspond to a ~ 0.3, Ra — 7 fm, and g"^^ ~ 2 GeV. Now that 
the model definition and the model parameters are given, we can calculate 
af^ and af^ numerically. Then, in the high-energy limit, it is a common 
assumption that the nucleus source is infinitesimally thin, i.e. 

Pt{x^,x~) = pt{x^)5{x~) , pp(a3_L,a;"^) = ^p(a;±)(5(x+) . (32) 

Then, the Wilson lines are replaced, respectively, by 

V\xi_) ^ e'^'^"'^^^) , W\xi_) ^ e'5^*'^(^^) (33) 

with the solution of the Poisson equation, 

- dlK^'\x^) = -p,{x^) , -9iA(P)(a;^) = -p^{x^) . (34) 

One should be very careful about this replacement because this is not an 
approximation on Eqs. (jlSp and ()2ip . Even though the longitudinal extent 
in the color source is infinitesimally thin, it should not be legitimate to drop 
the path ordering [22]. Therefore, we have to think that the numerical MV 
model is something distinct from the original MV model in the continuum 
variables. 

In the MV model the distribution of the color source is random and there 
is no correlation between different sites. Figure H] illustrates an example of 
the initial h.{x±) as a solution of the Poisson equation. Because the operator 
^/d\ involves spatial average, we see that the spatial distribution of A(a;^) 
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Fig. 5. Left: Time evolution of the chromo-electric and the chromo-magnetic fields. 
The subscripts, L and T, represent the longitudinal and the transverse fields, re- 
spectively. Right: The longitudinal and the transverse pressures as a function of 
time. 



is rather smooth even though the source p{x±) has random fluctuations. 
This smoothness is not physical, however, and the gauge fields are furiously 
fluctuating as shown in the right panel of Fig. HI It is worth noting that 
the color-flux tube picture as sketched in Fig. [3] is not the case in the MV 
model and the JIMWLK evolution is indispensable to take account of the 
flux tube structure. 

The physical observables are measured by taking an ensemble average 
of results with different initial p{x±ys. It is useful to compute not only the 
energy density pip but also other combinations of the energy-momentum 
tensor. In particular the following pressures are important in order to judge 
how anisotropic the system is; 

= i(T^- + Tyy) = (tr [El+Bl]) , (35) 

P, = (r'T^r^) = {tv[El + Bl - El - Bl]) , (36) 

where the longitudinal and transverse chromo-electric and chromo-magnetic 
fields are defined as 

El = g^'^g'T' , El = ^ (^gxagxa ^ gyagya^^ ^ (37) 

T32 77a 77a 732 ^ / 77a 77a I 77a 77a N {'JO\ 

~ -^12-^12 ) -"t ~ ^2 y-^vx-^vx vy vyJ ■ ^ ' 

The numerical results from the numerical Glasma simulation are pre- 
sented in Fig. [5j From this figure it is clear that there are only longitudinal 
fields El and Bl right at the collision (r = 0) as explained with Fig. [3l The 
transverse fields are developing as r increases, and eventually the longitu- 
dinal and the transverse fields approach each other at g^/xr > 1. This does 
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not mean isotropization, however, because there are two {x and y) compo- 
nents in the transverse direction. One can understand what is happening 
by taking a careful look at the pressures. The longitudinal pressure goes 
to zero for g^fiT > 1 and this means that particles move on the free stream 
along the longitudinal direction. Thus, the system remains far from thermal 
equilibrium. 

Here we encounter a perplexing situation. We know the the CGC should 
be a correct description of the initial dynamics of the heavy-ion collision. 
Even if the CGC cannot reach isotropization and thermal equilibrium, it 
should be a natural anticipation that the CGC can at least capture a correct 
tendency toward thermalization. This anticipation is not the case at all, 
however, as seen in Fig.[5j Then, is there anything missing in our discussions 
so far? 



Yes, there is. We have neglected fluctuations on top of boost-invariant 
CGC-background field and thus rj dependence at all. Such a treatment is 
not always justified. As a matter of fact, longitudinal structures in a finite 
extent by not using pt{x±^) 6{x~) but treating pt{x±,x~) correctly in the 
path ordering would give rise to r/-dependent fluctuations. Also, quantum 
fluctuations have rj dependence as well [H [23] . 

One might have thought that small disturbances in the longitudinal 
direction could make only a slight difference. But, the fact is that there 
is a tremendous difference between results with and without r^-dependent 
fluctuations. With random fluctuations in r/ space, the Glasma simulation 
would lead to significant decay from the CGC background fields with rj- 
independent zero-mode into ry-fluctuating non-zero modes [3]. An example 
is shown in Fig. [6l 

To obtain the results as presented in Fig. El only u = zb^o rnodes are 
disturbed. 



where is the longitudinal extent that we took as = 2 in the simulation. 
Once 6E^ is given, the transverse fluctuations, 5E^, are chosen in such a 
way to satisfy the Gauss law. Then, physical observables of our interest are 
Fourier transformed from rj space to i' space. It is now obvious that the 
CGC background fields at = are relatively larger than fluctuations by 
small A and the second dominant mode should sit at u = uq. In this article 
we limit ourselves to the simplest choice of uq = 1. 

Usually unstable modes grow up exponentially, but in the expanding 
geometry, the instability implies a slightly weaker growth according to ~ 



5. Glasma instability and the scaling spectrum 




(39) 
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Fig. 6. Left: Longitudinal pressure dX v = vq = 1 for various system sizes. Right: 
Spectrum of the longitudinal pressure at late time presented in v space. 
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V 

Fig. 7. Left: Instability as a function of the magnitude of the seed. Saturation 
behavior in the non-linear regime is seen slightly for A = 10^"^ and clearly for 
A = 10^^. Right: Corresponding spectra in the v space. 



exp[a-v/r]. The horizontal axis in the left panel of Fig. [6] is, thus, not the 
dimensionless time g^ixT itself, but \J g^\iT. We can surely confirm that the 
longitudinal pressure component at = fo = 1 increases linearly in Fig. [6] 
when plotted -with the logarithm of the pressure as a function of \J q^ilt. 

The detailed structure of the instability in the spectral pattern is inter- 
esting to see. The right panel of Fig. [6] is the spectrum corresponding to 
the simulation results shown in the left. 

From now on, let us consider the fate of the instability; it is simply 
impossible for the instability to last for ever. The saturation of instability 
growth can be observed in two ways. The first case is that the initial magni- 
tude A (appearing in Eq. ([5U|) ) is taken to be large enough to accommodate 
non- linear effects. The second is the large-time behavior - simply we wait 
until the unstable modes grow up considerably. 

As shown in the left panel in Fig.[7]the instability for A = 10~^ stops and 
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Fig. 8. Saturation of the instability at late time with indications to three distinct 
regions; the linear regime where the instability persist, the non-linear regime where 
the energy flow is expected, and the UV-cutofF regime where artificial results by 
the UV cutoff are unavoidable. 



the spectrum is flattened after all. In the saturated regime at \J g^[iT > 20 
we see that the pressure at u = i^q slightly decreases rather than increasing. 
This behavior can be interpreted as follows; As long as the non-zero modes 
are small (i.e. in the linear regime), the energy decay from the dominant 
CGC background at v = makes non-zero modes enhanced exponentially. 
The injected energy is much bigger than the escaped energy toward higher-i/ 
modes in the linear regime. This balance changes gradually with increasing 
amplitude of unstable modes, and eventually a steady energy flow is ex- 
pected when non-Abelian nature of non-zero modes becomes substantially 
large (i.e. in the non- linear regime). At even larger time, as hinted from 
Fig- 13 the spectrum is flattened and the UV-cutoff effects at large v should 
be appreciable. Then, an intriguing question is; what is going on in the 
non-linear regime before the UV-cutoff effects contaminate the simulation? 

To address this question, we shall take the latter case, namely the long- 
time run of the simulation with small A for better numerical stability. The 
qualitative features in the results for A = 10~^ shown in Fig. [8] are just the 
same as the top curve in Fig. [7] that represents the results for A = 10~^. 
We are interested in the energy spectrum in the non-linear regime that can 
be immediately identified on Fig. [8l 

Figure [9] is the energy spectrum (not the Fourier-transformed longitu- 
dinal pressure but the energy contribution from respective ly modes). That 
is, what is shown in Fig. [9] is given as 

e^{u) = {tT[£'"'{-iy)£'^''{u) + T~^£"'{-u)£'^{u)] > . (40) 

We see that there is a clear tendency to approach a scaling form in the 
non-linear regime (as in the right panel of Fig. E]). With some rescaling 
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Fig. 9. Left: Energy spectrum approaching the power-law scahng oc v Right: 
Breakdown of the Kolmogorov-type scahng in the UV-cutoff regime. 



we find that this power is exactly consistent with the Kolmogorov value, 
—5/3. Generally speaking, in non-expanding systems, it is not a surprise 
that the Kolmogorov power emerges because this value can be guessed by 
dimensional analysis. Such a dimensional argument may work even for the 
expanding system. In the Bjorken coordinates r] is dimensionless, but the 
physical scale in the longitudinal direction is to be interpreted as rr/, and 
thus the corresponding momentum should be u/t. Then, this quantity 
gives a dimension of length (with c multiplied appropriately). Because of 
an expected energy flow in u space, its rate '0 is also a relevant quantity. 
Then, the characteristic length and time scales of the system are expressed 
by two quantities with the following dimensions; 

[v/t]=1-^ , ['il;] = fr^. (41) 

The energy spectrum, on the other hand, has the dimension, 

[V^T'eA^)\=lH-\ (42) 

that is reproduced exactly by a unique combination of (z^/r)~^/'^(^/^)^/^. 
Therefore, it is concluded that ^{v/t) should exhibit the power-law scal- 
ing in terms oiu/r whose power is identical to the Kolmogorov's value, —5/3 
as long as the system stays in the non-linear regime. From the left panel 
of Fig. [9l we can see only the scaling region or the so-called inertial region 
realized. The dissipative region at high u is not found, probably because of 
the UV-cutoff effects. 

This nice scaling is lost at further later time. We can understand from 
the right panel of Fig. [S] that the energy flow is stuck at the UV edge and the 
energy spectrum is artificially pushed up by the UV-cutoff contamination. 
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6. Outlooks 

It was a surprise that the Kolmogorov's scahng law could be realized 
in the Glasma simulation in the expanding geometry. The dimensional 
argument is not so strong to constrain the shape of the energy spectrum 
uniquely. It should be an important test whether the Kolmogorov's —5/3 
behavior can be confirmed or not in other simulations of the pure Yang- 
Mills theory. It has been established that the non-Abelian plasma generally 
has an instability associated with anisotropy that grows up until the non- 
Abelianization occurs |24] . Similar phenomena of instability tamed by non- 
Abelian interactions are found in many other examples. Then, presumably, 
there must appear the power-law scaling in the region around the non- 
Abelianization. 

The turbulent decay and the associated Kolmogorov's power-law are 
interesting discoveries from the long-run simulation of the Glasma. But, 
it cannot answer anything about the realistic thermalization mechanism 
from the Glasma. The turbulence is certainly a tendency into thermalized 
matter, but the energy flow is a slow process and the relevant time scale, 
g^HT ~ 2000, is far outside of the validity region of the CGC description. 
There must be still something missing that can accelerate the thermalization 
speed. 

If this missing piece were finally set in, the transient Glasma would 
provide us with the initial input for the hydrodynamic equations. Even in 
this case, the analysis we have seen should be useful. The energy spectrum 
in V space should carry important information. Then, the power could agree 
with or deviate from —5/3 as suggested in strong-coupling studies pSj |26]. 
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